#work out whether SNPs dn or ds
#reworked version to use Biostrings not GeneR

#wrapper for a per contig analysis
codingSNPref <- function(AllSNPref,dir.fasta="/Paterson/Datafiles/grouse/gContigs",dir.embl="/Paterson/Datafiles/grouse/embl_exon",seq.name=NULL,file.name=NULL){
	#AllSNPref is an object (dataframe) from baseCoverage
	#expects embl files in dir.embl, perhaps created by add.CDS
	
	if(is.null(AllSNPref)) return('null.input')
	if(identical(AllSNPref,0)) return('zero.input')
	
	
	if(is.null(seq.name)) seq.name <- AllSNPref[1,"cDNA"]
	if(is.null(file.name)) {
		fasta.file.name<-paste(dir.fasta,"/",seq.name,"_matched.fasta",sep="")
		embl.file.name<-paste(dir.embl,"/",seq.name,"_exon.embl",sep="")
		}else{
		fasta.file.name <- file.name[1]
		embl.file.name <- file.name[2]
		}
	
	if(!file.exists(fasta.file.name)){ #would this ever be called?
		warning(fasta.file.name, " does not exist... returning\n")
		AllSNPref$cds <- FALSE
		AllSNPref$dn <- FALSE 
		return(AllSNPref)
		}
	if(!file.exists(embl.file.name)){
		warning(embl.file.name, " does not exist... returning\n")
		AllSNPref$cds <- FALSE
		AllSNPref$dn <- FALSE 
		return(AllSNPref)
		}

	require(Biostrings)
	
	#for some reason won't read straight into buffer
	ref.seq.set <- read.DNAStringSet(fasta.file.name)
	ref.seq <- DNAString(ref.seq.set[[1]])
	
	#extract ferature table
	tmp.features <- system(paste("grep \'FT\'",embl.file.name,"|sed \'s/FT[ \t]//\'"),intern=TRUE)
	cds.lines <- grep('^[ ]*CDS ',tmp.features)
	product.lines <- grep('/product',tmp.features)
	cds.string <- character(length=length(cds.lines))
	for(cdsi in 1:length(cds.string)){
		grep.cds <- paste(tmp.features[cds.lines[cdsi]:(product.lines[cdsi]-1)],collapse="")
		
		}
	#extract coding seqs
	grep.cds <- sub('^[ ]*CDS ','',grep.cds)
	grep.cds <- gsub('[ ]*','',grep.cds)

	#grep.cds <- system(paste("grep \'CDS\'",embl.file.name,"|sed \'s/FT[ \t]*CDS[ \t]*//\'"),intern=TRUE)
	grep.cds <- sub("join\\(","",grep.cds)
	grep.cds <- sub(")","",grep.cds)


	AllSNPref$cds <- FALSE
	AllSNPref$dn <- FALSE
	AllSNPref[,c('cds','dn')] <- t(sapply(X=AllSNPref[,"id"],FUN= catagoriseSNPs,grep.cds=grep.cds,ref.seq=ref.seq))
	AllSNPref
	}



#seq.name <- "cDNA_260-1"
#dir.embl <- "/Paterson/Datafiles/grouse/embl_exon"

#read sequence into buffer




#sequence 0 for original seq, seq 1 for mutated seq, seq 2 for temporary seq
#translateCDS <- function(cds.str,seqno=0){
	#expects cds.str in for 1..10,20..30 etc
	#and for a sequence in buffer seqno
#	tmp.from <- as.numeric(sapply(X=strsplit(cds.str,","),FUN=function(X){sub('\\.\\.[0-9]*$','',X)}))
#	tmp.to <- as.numeric(sapply(X=strsplit(cds.str,","),FUN=function(X){sub('^[0-9]*\\.\\.','',X)}))
#	assemble(seqno=seqno,from=tmp.from,to=tmp.to,strand=0,destSeqno=2)
#	translate(2)
#	}
	
mutateSeq <- function(snp.id,cds.str,ref.seq){
	#expects a snp.id, ie id from baseCoverage
	#cat(snp.id,"  ")
	snp.id <- sub("^cDNA_","cDNA--",snp.id)
	snp.id <- sub("^ssh_","cDNA--",snp.id)
	tmp.from <- as.numeric(sapply(X=strsplit(cds.str,","),FUN=function(X){sub('\\.\\.[0-9]*$','',X)}))
	tmp.to <- as.numeric(sapply(X=strsplit(cds.str,","),FUN=function(X){sub('^[0-9]*\\.\\.','',X)}))
	tmp.id <- strsplit(snp.id,"_")[[1]]
	#if(substr(tmp.id[1],2,6)=="ontig"){
	#	mut.start <- as.numeric(tmp.id[2])
	#	}else{
	mut.start <- as.numeric(tmp.id[2])
	#	}
	
	
	mut.cds <- FALSE #test whether in coding region
	mut.dn <- FALSE
	if(any(mut.start>tmp.from&mut.start<tmp.to)){
		mut.cds <- TRUE
		}else{
		#if the mutation isn't in the coding region don't bother
		#working out anything else
		return(list(mut.cds=mut.cds,mut.dn=mut.dn))
		}

	mut.old <- sub("-","",tmp.id[5])
	mut.new <- sub("-","",tmp.id[6])
	if(nchar(mut.old)!=nchar(mut.new)){
		#if there's an indel in the coding sequence assume non-synonymous and return
		mut.dn <- TRUE
		return(list(mut.cds=mut.cds,mut.dn=mut.dn))
		}

	#should now be carrying on only if a straight swap of nucleotide(s)
	mut.len <- nchar(mut.old)
	mut.stop <- mut.start + mut.len -1

	#will probably fail for large sequences
	cds.old <- character(length=length(tmp.to))
	
	#must be quicker way, but generally only 1 to 3 exons
	for(cds in 1:length(tmp.to)) cds.old[cds] <- substr(as.character(ref.seq),start=tmp.from[cds],stop=tmp.to[cds])
	cds.stringset <- DNAStringSet(cds.old)
	
	#can't get translate to work on a stringset for some reason
	old.prot <- sapply(cds.stringset,function(X){
		translate(maskMotif(X,'N'))
		})

	#comes from Biostrings, only works for equal length swaps
	#could presumably save some effort by using "-" as a letter
	new.seq <- replaceLetterAt(ref.seq,at=seq(mut.start,mut.start+mut.len-1),letter=mut.new)
	
	cds.new <- character(length=length(tmp.to))

	for(cds in 1:length(tmp.to)) cds.new[cds] <- substr(as.character(new.seq),start=tmp.from[cds],stop=tmp.to[cds])
	cds.stringset <- DNAStringSet(cds.new)
	new.prot <- sapply(cds.stringset,function(X){
		translate(maskMotif(X,'N'))
		})
	
	if(!identical(sapply(old.prot,as.character),sapply(new.prot,as.character))) mut.dn <- TRUE
	
	list(mut.cds=mut.cds,mut.dn=mut.dn)
	}


catagoriseSNPs <- function(snp.id,grep.cds,ref.seq){
	tmp.mut.list <- vector(mode="list",length=length(grep.cds))
	for(i in 1:length(grep.cds)){
		tmp.mut.list[[i]] <- mutateSeq(snp.id,grep.cds[i],ref.seq)
		}
	c(any(sapply(tmp.mut.list,function(X){X$mut.cds})),
		any(sapply(tmp.mut.list,function(X){X$mut.dn})))

	}




